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ABSTRACT 



An extensive review of the literature revealed that many 
experimental studies have been conducted in heat-transfer wind 
tunnels simulating the cooling passages in turbine blades. 
However, very few numerical studies have been performed. 
Phoenics, a computational fluid dynamics computer program, 
produced results for several duct configurations and calculat- 
ed the heat transfer characteristics of each. The configura- 
tions investigated included a straight, square duct and a 
rectangular duct with turbulators (a form of turbulence 
promoters) present. The parameters varied included Reynolds 
numbers, turbulence intensity and grid geometry. Results for 
the turbulated duct indicated highly distinct and repeatable 
flow patterns developed over a wide range of values. The 
variation of the inlet turbulence intensity had little impact 
on the kinetic energy. The results proved to be highly grid 
dependent, which greatly impacted the correlation between 
experimentally and numerically produced data, for the same 
configurations, under similar operating parameters. 
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I . INTRODUCTION 



A. BACKGROUND TO TURBINE INTERNAL COOLING 

All modern gas turbine engines share three basic 
components; the compressor, the combustion chamber and the 
turbine. From the compressed mainstream flow out of the 
combustion chamber the turbine extracts the power to operate 
the compressor and to deliver the gases to either a power 
turbine for shaft work or to a nozzle for thrust . During 
this process losses occur in both the compressor and the 
turbine sections. In the compressor they serve to increase 
the power absorbed while in the turbine they act to decrease 
the power output. To overcome these losses, more air and 
therefore fuel flow are required for a required power level. 
The compressor loss increases the temperature of the flow out 
of the compressor and into the combustor and decreases the 
amount of fuel which can be added for a given turbine inlet 
temperature. Based on cycle analysis, the higher the 
allowable temperature into the turbine the higher the power 
output of the unit. The turbine inlet temperature constrains 
the design of all gas turbine engines. 

Thus for a given engine, the amount of fuel which the fuel 
control unit (FCU) can add is limited by the ability of the 
materials to withstand thermal loads. Failures can result 
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from effects such as oxidation, corrosion and erosion, as well 
as creep (due to a combination of thermal and centrifugal 
loads) and fatigue, which ensues from both thermal and cyclic 
loads. The combination of these effects determine the 

operating life of the turbine. Two solutions to this problem 
have emerged. One is the improvement of materials for the 
construction of turbine components, the other is the use of 
compressor bleed air to cool the turbine components. 

Turbine blade cooling will be discussed in the following 
order: external cooling air over turbines, impingement 

cooling, film cooling, transpirational cooling and with the 
final focus on internal cooling using turbulators . 
Experiments in turbine cooling began in the 1940 's with the 
first units going into operation around 1960. Initial engines 
bled 1.5 to 2 % of the compressed air from the compressor to 
cool the turbine 200-300° C [Ref 1] . Modern aircraft engines 
can use as much as 20% of the compressor discharge flow [Ref 
2] . This cooling air externally bathes the turbine blade with 
initial rotors receiving their cooling air directly from the 
compressor discharge. Two shortcomings associated with this 
type of cooling include the turbine's reduced ability to drive 
the compressor because of the lower temperature of the flow as 
it enters the turbine and the mixing of the cooling and the 
mainstream flow that results in increased aerodynamic losses. 
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Rather than just providing external cooling air, 
individual cooling of the airfoil itself is another attempt to 
elevate turbine inlet temperatures. Impingement cooling at 
the leading edge and at midchord consists of placing a sheet 
metal insert into the hollow turbine blade with appropriately 
placed venting holes . Coolant flows through the root of the 
blade into the insert and out the impingement holes. This 
convective type cooling directs air at the hottest spots on 
the blade and creates airflow inside the blade. Often the 
most difficult part of the design analysis involves 
determining the exact location of the stagnation point or 
hottest spot on the blade. After the cooling process occurs 
the blade vents the coolant either at the tip or at the 
trailing edge, which then joins the main flow. 

Film cooling also uses air from the compressor routed into 
the root of the blades. In film cooling this air travels 
through the blade and departs through openings designed to 
produce a blanket of cooling air to protect the blade against 
the mainstream gases . The analysis on the design of these 
holes focuses on number, location, size, shape and orientation 
angle. With the average cooling effectiveness (ACE) defined 
as ; 
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where T tg is the mainstream gas total temperature, T m is the 
average surface temperature of the airfoil and T c is the 
coolant supply total temperature, engineers usually specify 
film cooling for an ACE of 0.6. Film cooling is very 
effective when combined with other forms of internal 
convective cooling, but the price in total pressure losses for 
putting such a quantity of cooling air into the hotter main 
gas flow is considerable. 

Another form of cooling is transpirational or aspiration 
cooling. This technique uses a porous material to construct 
the airfoil and runs cooling airflow into and through it. 
This combines both convective and film cooling, but with some 
obvious problems associated with the strength of the airfoil 
construction. Designers using a woven wire material achieved 
some success, however difficulties in the shaping of the 
airfoils may preclude this type of construction. 

The final type of cooling discussed is internal cooling 
passages with coolant run through the blade and out at the tip 
and trailing edge of the blade. Only convective cooling is 
used and chosen when the design calls for an ACE of less than 
0.5. Internal designers choose convective cooling because of 
a limited supply pressure (a higher ACE would require a higher 
pressure) , and a higher effectiveness level using just 
convective cooling could lead to a higher temperature gradient 
and thus aggravate the very thermal stress problem which 
cooling tries to avoid. 
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B . EXPERIMENTAL INVESTIGATIONS 



Investigations into internal passages for convective heat 
transfer fall into three sub-areas; straight passages, 
serpentine passages with one or more 180° bends and single and 
multi-pass channels undergoing rotation. The investigations 
often include analysis of different flow properties, channel 
shapes, aspect ratios and wall roughness, which encompasses 
the shape, location, number and orientation of turbulence 
inducing structures. The continual refinement of the 
operating characteristics of turbine engines has made this 
last characteristic increasingly important as many of the 
other parameters are fixed. 

As early as the 1960 ' s research centered on the effects of 
protrusions into the flow from the walls of the channel. 
Wilkie [Ref 3] was one of the first to publish results of 
experimentation on the effects of rib roughened ducts on heat 
transfer. The interest in cooling operations in nuclear 
reactor fuel elements resulted in the findings he collected 
and published. Through variations in duct equivalent diameter 
(d) , rib height (e) and rib pitch or spacing (p) he found 
correlations between the Stanton number, the Reynolds number, 
the p/e ratio and the e/d ratio. The Stanton number is 
defined as h/ (upC p ) where h is the heat transfer coefficient, 
u is the velocity in the free stream at the outer edge of the 
boundary layer, p is the fluid density and C p is the specific 
heat at constant pressure. The results showed Stanton number 
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varied less with Reynolds number for small values of e/d 
(c.012), and behaved much as on a smooth surface at larger 
values of e/d (>.012) . In addition he found the value of the 
p/e ratio giving maximum Stanton numbers varied with e/d and 
Reynolds number but tended to approximately the value six at 
large e/d irrespective of Reynolds number. Norris [Ref 4] 
also detailed the effect of these " turbulators " in his 
investigation. After testing various shapes, he concluded 
that although the local heat-transfer coefficient varied in a 
very non-uniform manner in the stream-wise direction, there 
was a value of the pitch-to-width ratio for the square 
turbulator which gave an optimum value for the average heat 
transfer as well as friction factor. Increasing the ratio 
above this level served only to decrease these values. Also, 
he developed a relationship between the friction factor of his 
test apparatus and the friction factor of a smooth duct, and 
found an optimum value existed for this ratio. Increases 
above this value did not give increased heat transfer. 

Using the Nusselt number as an indicator of heat transfer, 
Han, et . al . [Ref 5] varied their rib angle (a) and spacing to 
determine their effects. With a rib-to-height ratio of ten, 
they found the Nusselt number to be twice as high with the 
turbulators placed 90° to the flow as for the case without 
turbulators. Introducing an a of 60° increased the Nusselt 
number another 5% while an a of 45° delivered an additional 
25% over the 90° orientation. Lau, et . al . [Ref 6] built upon 
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this concept when they examined fully developed flow in a 
square channel with turbulence promoters on two opposed walls, 
using turbulators configured with a's of 45°, 60° and 90°, 
lengths running from wall to wall, and shorter discrete 
turbulators. Both configurations were tested, crossed or 
parallel with respect to each other. They found parallel 
discrete ribs with a's of 60° provided the highest ribbed wall 
heat transfer, while the same style of ribs canted at 30° 
provided the lowest level of pressure drop. In general they 
recommended using parallel, angled discrete ribs for internal 
cooling passages as they provided higher ribbed wall heat 
transfer, lower smooth wall heat transfer and lower channel 
pressure drop. The optimum configuration for these discrete 
ribs was a continually staggered row pattern of three ribs 
followed by a row of two ribs. 

Taslim and Spring [Ref 7] centered their straight square 
channel experiments not so much on the placement of the 
turbulators, but rather the geometry of them. They determined 
that an optimum turbulator spacing existed for every given 
turbulator blockage ratio and aspect ratio. The turbulator 
blockage ratio was defined as the ratio of the turbulator 
height to the hydraulic diameter of a non-turbulated section. 
They also found that the Nusselt number decreased as the 
turbulator aspect ratio decreased. This became even more 
critical as the spacing dropped below the optimum value. They 
also developed an alternatively shaped turbulator which not 
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only provided the same heat transfer performance as a larger 
turbulator, but also enhanced the castability and, thus, 
reduced the probability of problems in manufacturing. To do 
this they developed a turbulator with a narrow tip width to 
enhance flow while at the same time widening the base to ease 
the problems of castability. The resulting shape is 
trapezoidal with slightly rounded corners at the tip. 

Prior to Boyle's [Ref 8] report, most data for turbulators 
came from experiments using long passages with measurements 
done in regions where results did not change with distance. 
Boyle initiated his report with an investigation of a 
serpentine channel without any turbulators. He researched the 
effect of bends only on heat transfer coefficients. He found 
that when varying the outside edge of the turn from square to 
a round geometry he did not alter the heat transfer rate. He 
established that the rate varied directly proportional to the 
Reynolds number and the highest rate occurred at the beginning 
of the second leg after the bend. The addition of turbulators 
markedly increased the heat transfer rate. He also found an 
optimum height which required the turbulator to extend into 
the boundary layer buffer region. 

Chandra, et . al. [Ref 9] did additional work on a multi- 
pass serpentine channel investigating the effect of the 
turbulator angles on the flow. Using a ratio of average 
Sherwood numbers normalized to Sherwood numbers for fully 
developed turbulent flow, they equated this value to the 
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ratios of average Nusselt numbers normalized with respect to 
Nusselt numbers for fully developed turbulent flow. Here 
they define the Sherwood number (Sh) as: 

Sh= httPh (2) 
VS C 

where h^ is the local mass transfer coefficient or the local 
mass transfer rate per unit area divided by the difference of 
the coolant density at the wall and freestream; D h is the 
hydraulic diameter, v is the kinematic viscosity and S c is the 
Schmidt number for the coolant . The Schmidt number is the 
ratio of the kinematic viscosity to the mass diffusivity. 
Using the above equation they were able to make a heat to mass 
transfer analogy. For the straight channel they found 
periodic Sherwood numbers on the ribbed walls in the direction 
of flow, while a's of 60° and 45° produced higher Sherwood 
numbers on the outer walls than on the inner ones. 
Turbulators set perpendicular to the flow produced little 
spanwise variation of heat transfer, and the slight spanwise 
variations they did find diminished somewhat after the 180° 
turns. In general they found that with respect to the area of 
the turn the overall average Sherwood number ratio was 
independent of the rib angle but tended to decrease slightly 
as the Reynolds number increased. 

Using a three-pass rib-roughened channel Han and Zhang 
[Ref 10] established the local pressure drop distribution. 
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With these data, they were able to determine the local coolant 
flow rate distribution. Their results indicated that the 
pressure dropped linearly for the first leg to the turn. Just 
prior to the turn, a slight rise followed by a sharp drop 
occurred into, through and just after the bend. Flow 
recirculation caused a recovery of pressure during the second 
straight channel, with a similar pressure drop occurring at 
the second turn. They also found the turbulators' biggest 
impact to be evident in the first straight channel prior to 
the turn. The 180° turns tended to dilute their impact on 
subsequent legs . 

Chyu [Ref 11] did similar testing with both two- and 
three-pass channels employing 180° turns. The emphasis of 
these experiments was on surface heat transfer coefficient 
versus pressure distribution. With the two-pass channel, he 
found that while the pressure loss attributable to the turn 
was significant, it was independent of the Reynolds number. 
His findings of flow separation, reattachment, impingement and 
secondary flows caused by the 180° turn are consistent with 
the findings of Han and Zhang. He attributed the heat 
transfer levels observed to be due to these secondary flow 
effects. Also, even though he found the heat transfer rates 
of both turns to be approximately equal, the total heat 
transfer of a three pass channel was about 15% greater than 
for a two pass channel. 
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With a three-pass square channel, with two 180° turns, Han 
and Zhang [Ref 12] investigated the impact of angling 
turbulators 60° on opposite walls. They also opted for mass 
transfer to be the measured experiment parameter and converted 
to Sherwood numbers to correlate the heat transfer. 
Configuring the ribs both crossed and parallel to the opposite 
wall, they found a's of 60° produced greater mass transfer 
coefficients than a's of 90°, while parallel ribs generated 
greater coefficients than crossed ones. Starting with a 
smooth channel, they established the value of the second leg's 
mass transfer to be 60% greater than the first leg. Since the 
third leg closely approximated the second, they surmised that 
the 180° turn caused the difference. Finding, in general, the 
60° parallel orientation provided the highest mass transfer 
coefficients, this study showed rib angle orientation plus 
sharp turns in the coolant path can increase or decrease the 
mass transfer coefficient. In addition, by correlating 
previous data with the results of this experiment, they 
validated the application of straight ribbed channel data to 
the first pass in multi-pass ribbed cooling channels. 

The final area of investigation associated with internal 
convective cooling of turbine blades relates to rotating 
elements. Clifford [Ref 13] used a rotating, three pass, 
smooth channel with airflow directed through it. Measuring 
the heat transfer along the leading and trailing walls with 
respect to rotation, he found rotation served to both aid and 
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inhibit heat transfer with the reductions approaching 30% on 
the suction side of the final leg and the enhancement closer 
to 35% on the pressure side of the initial leg. Averaging 
his data showed the effects of rotation reduced the heat 
transfer to 90% of the statically measured values. 

Morris and Harasgama [Ref 14] performed similar tests 
using axial airflow directed both radially inward and outward. 
Their results showed outwardly directed airflow promoted 
greater cooling on the trailing (pressure) surface vice the 
leading (suction) surface. Inwardly directed airflow, 
however, increased heat transfer on the leading surface. This 
heat transfer was rotational Reynolds number dependent as it 
increased slightly as the number was increased. Like 
Clifford, they also discovered instances of reduced heat 
transfer when compared with the stationary tube condition. 
They also found rotational buoyancy tended to suppress heat 
transfer with radially outward flow. With respect to 
stationary experiments they found no areas of impaired heat 
transfer . 

The findings of Kheshi and Scriven [Ref 15] and Harasgama 
and Morris [Ref 16] were consistent. They found that as the 
Rossby number (which is the ratio of the mean through-flow 
velocity of fluid along the duct to a peripheral velocity 
formed by the product of the angular velocity and the 
hydraulic diameter (Qd/V) ) increased, the axial velocity on 
one side decreased with respect to the other side and a cross 
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channel flowfield developed which was perpendicular to the 
main flow. The steady state for this condition showed that a 
two-vortex pattern developed. Increasing the Rossby number 
further generated a four-vortex flow pattern in the steady 
state. This four-vortex pattern required a shear-free flow to 
develop on both the top and bottom walls. The experiment 
suggested the flow pattern was a function of the aspect ratio 
of the channel, the Rossby number and the Ekman number. The 
Ekman number is defined as the ratio of the Rossby number to 
the Reynolds number [Ref 17] . Harasgama and Morris confirmed 
the effect of heat transfer on both the leading and trailing 
edges with either outward or inward flow. In addition they 
put a value of 10-20% on the decrease of the heat transfer on 
the leading side during outward flow. They attributed the 
increase in trailing side heat transfer to the strengthening 
of secondary flow as the average Nusselt number increased. 

Hajek et . al . [Ref 18] investigated the impact of the rate 
of rotation and the difference in the temperature between the 
coolant and the wall (AT) . By increasing AT the heat transfer 
either remained constant or increased. This increase was 
greater for walls with turbulators, although the percentage 
was not. This was a result of the higher heat transfer to 
walls with turbulators. When these increases, which could be 
larger, were compared to their larger initial values, smaller 
percentage increased were found. Increasing the Rossby number 
induced higher heat transfer on the trailing surface while at 
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the same time decreasing the leading surface value. They also 
developed the buoyant parameter (BP) such: 



BP=(^E) (f^i) (£^) 

p V V 



(3) 



where Ap is the pressure differential, p is the density of the 
coolant, Q is the angular velocity, W is the eccentricity of 
the mid cross-sectional plane of the duct, V is the mean 
coolant velocity along the duct and D h is the hydraulic 
diameter of the duct. By increasing this value, a higher heat 
transfer occurred on the trailing surface, while the leading 
surface initially experienced a decrease with a subsequent 
increase when BP approached the 0.2 to 0.4 range. In general 
they found that the locus of maximum heat transfer ratios for 
each Rossby number was a function of BP. 

Taslim, Rahman and Spring [Ref 19] also confirmed that an 
increase of heat transfer occurred for both the rotating and 
stationary cases in a single pass channel when staggered rows 
of turbulators were added on opposite walls. The main thrust 
of this experiment was to explore the effect of Coriolis 
forces on internal heat transfer. In addition, they 
investigated the effect of varying the blockage ratio, e/D h , 
where e denotes the turbulator height and D h designates the 
hydraulic diameter of the passage. Their experiments found a 
blockage ratio of 0.1333 produced a maximum heat transfer 
coefficient increase of 45% and a ratio of .333 generated a 
minimum, a 6% decrease. Guidez [Ref 20] used a rotating 



14 



channel to correlate Reynolds number and Nusselt number to 
heat transfer coefficient. His results indicated that the 
pressure (or trailing) side of the channel experienced 
slightly increased, heat transfer as the Reynolds number 
increased while the suction side heat transfer was weaker. He 
determined the turbulence at the higher Reynolds number became 
larger than the turbulence induced by the rotation. His ratio 
of rotating Nusselt number to the stationary Nusselt number 
increased from 1.0 to 1.75 as the channel accelerated from 0 
to 600 radians/sec on the pressure side. The ratio on the 
suction side decreased from 1.0 to 0.7 and reached a steady 
state value at approximately 175 radians/sec. 

Wagner et . al . [Ref 21] performed an in-depth experiment 
to determine the effects that forced convection, Coriolis, 
buoyancy and flow direction had on heat transfer. For outward 
flow, changes in density ratio and Rossby number significantly 
affected the heat transfer. Altering the buoyancy parameter 
for this flow caused large changes in the heat transfer for 
the high pressure side, while changes made during inward flow 
produced minimal changes in the heat transfer. Density ratio 
and Rossby number also contributed minimally during inward 
flow. The suction side of the channel proved relatively 
immune to changes in heat transfer caused by modifying the 
density ratio, Rossby number or the flow direction. The heat 
transfer was more a function of the buoyancy than any other 
parameter . 
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Morris and Hanami-Nasr [Ref 22] found similar results in 
the actions of outwardly directed flow. They found enhanced 
heat transfer on the pressure wall resulted from Coriolis 
induced secondary flow and a reduction on the suction face. 
They determined that turbine blades designed from heat 
transfer measurements made on stationary models could 
experience failures attributable to the inability to predict 
local hot spots . 

Investigating the effects of channel and turbulator 
geometry on internal heat transfer, Taslim et . al. [Ref 23 and 
24] investigated a spanwise rotating channel. In both of 
these experiments, two walls of a channel contained 
turbulators arranged at an angle of attack to the direction of 
flow (a) of 45° in a criss-cross arrangement. In both 
experiments, blockage ratios, aspect ratios and Reynolds 
numbers were varied in an attempt to determine their effect on 
heat transfer coefficients. Using the Rotation number, which 
is defined as; 



Ro= 



OD h 

u m 



1 

’RoFsBy~JjumEer 



(4) 



where Q. is the angular velocity, D h is the hydraulic diameter 
and U m denotes the mean velocity of air, the results of both 
investigations showed rotational effects to be more pronounced 
when the test sections possessed higher aspect and lower 
blockage ratios. Specifically, a steady increase in heat 
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transfer coefficient observed on the trailing side as the 
Rotation number increased while the heat transfer coefficient 
on the leading side showed a steady decrease. 

C . COMPUTATIONAL INVESTIGATIONS 

Taylor et . al . [Ref 25], Moore and Moore [Ref 26], Medwell 
et . al . [Ref 27], Patankar and Prakash [Ref 28], Rowley and 
Patankar [Ref 29] and Patankar and Kelkar [Ref 30] all 
reported results relevant to the next phase of the present 
investigation: a computational fluid dynamic analysis of 
cooling airflow inside a turbine blade. All used a finite 
volume analysis which produced good correlation to either 
problems they ran concurrently or to rudimentary experiments 
performed previously. In addition, by including the conduction 
in the solid boundary and avoiding upwind techniques, they 
concluded that this procedure will prove to be exceedingly 
successful with more complex shapes in the future. 

D . PURPOSE 

Numerous experiments have been performed to investigate 
the flowfield associated with an internally cooled turbine 
blade. There have been considerably fewer attempts to analyze 
this flow using numerical techniques. The purpose of the 
present analysis is to use Phoenics, a PC-based computational 
fluid dynamics program, to assist in predicting and 
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understanding the flowfield based upon the experimentations of 
Taslim and Spring [Ref 7] . 
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II. THEORY 



A. GOVERNING EQUATIONS 

For a Newtonian fluid with no body forces, the law of 
Conservation of Mass becomes the continuity equation, Newton's 
Second Law yields the momentum equation and the First Law of 
Thermodynamics produces the energy equation. To correctly and 
fully describe fluid flow it is often necessary to include 
additional equations designed to model turbulence in the 
analysis. Finally additional equations are necessary to 
characterize the relationship between various fluid 
properties. An example of such an equation is the equation of 
state . 

1. The Continuity Equation 

When considering fluid flow through a fixed control 
volume, the continuity equation is: 

|£+v-(pv)=0 (5) 

ot 

where p is the fluid density and V with the overhead arrow is 
the total velocity vector and u and v are velocity components 
in the x and y directions for two dimensional flow. Using the 
substantial derivative 
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( 6 ) 



D 0 _ 0 () 

r>t dt 



+^-V() 



equation 5 becomes 



^£+p (V-v) =o 



(7) 



For incompressible flow 



Df> = 



Dt 



=0 



( 8 ) 



and equation 7 for incompressible two-dimensional flow becomes 

(9) 






dx dy 



2. The Momentum Equation 

The momentum equation for two-dimensional flow and in 
conservative form is [Ref 31] 



0pu + 


dp u 2 + 


dpuv 
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dx 
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=0 



( 10 ) 



( 11 ) 



Here the viscous effects of the flow are described using the 
shear stress terms T xx/ T xy and T yy where 
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( 12 ) 



2 du dv , 



^ * dy dx 



(13) 



| M 2|r-|^; 

3 dy dx 



(14) 



Now introducing Reynolds averaging, equations 10 and 11 
rewritten for incompressible flow assume the following forms 



|E + S,x ^*3(7 pPiP, (15) 

dt dx dy ox dx ** dy ** 



dp dp u v + dp v 2 __ dp f _d_ 
d t dx dy dy dx 



(^"P u 



V 7 



)+ -^^yy-<>v Tp) 



(16) 



where the reduced forms 



r »=“ (2 -|> 



(17) 



-xy 



^(|y + |T) 

dy dx 



(18) 



Vl“2^> < 19 > 

apply. In equations 10 through 19, u and v are the velocity 
components in the x and y direction respectively , p is the 
fluid density, p is the fluid pressure and [l is the fluid 
viscosity. All terms with an over bar indicate a time 
averaged quantity, while terms with a primed superscript 
indicate fluctuation quantities. 
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3 . The Energy Equation 

Applying the First Law of Thermodynamics to a fluid 
passing through a control volume yields 

pCL— =^P+V(icVT) +4> (20) 

H p Dt Dt 

where the total energy is comprised only of internal and 
kinetic energy. The form of equation 20 for heat transfer in 
incompressible two-dimensional flows is written in terms of 
static temperature as 



pC 4f tpC hl (u n * pC 4 |vT) = 



p£ + uPfc + vM + ( 21 ) 

dt dx dy dx dy 

4 (*f-P c p TFF> + <*§• -pc p Vf>) 

where the dissipation function, <I> overbar, is 



$ =— ( ) + — ( Pe + pl) +— ( PX) 

dx 1 ^ dy dx ' yy K dy’ 



( 22 ) 



and the fluctuating dissipation function has been ignored. 
[Ref 32] 

4 . Turbulence Modeling 

The unsteady Navier-Stokes equations are believed to 
describe turbulent fluid flow. Currently the primary 
technique employed in computational fluid dynamic predictions 
of fluid flow and heat transfer is to use the time-averaged 
Navier-Stokes equations. These take the randomly generated 
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flow variables and decompose them into their time averages and 
their instantaneous fluctuations. The time averages in this 
report are indicated with the overbar symbol while the 
fluctuations are given the apostrophe superscript. These new 
terms are interpreted as apparent stress gradients and heat- 
flux quantities associated with turbulent motion. These 
quantities can then be related to the mean flow variables 
through turbulence models. 



the two equation model, commonly referred to as the k-e model, 
discussed by Launder and Spalding [Ref 33] . The value for the 
effective turbulent viscosity, |a T , is 



where equals 0.09, k is the kinetic energy of turbulence 
and e is the turbulence dissipation rate. They derived two 
transport equations to define their turbulence model. The 
turbulence kinetic energy transport equation is commonly 
modeled as 



where Pr k is the Prandtl number for the turbulent kinetic 
energy and is assumed to be 1.0. The other transport equation 
necessary to complete this turbulence characterization is the 



For the present analysis, the turbulence model used is 




( 23 ) 



p|£ +p y.Vk=V( J^Vk) $-e 

dt Pr v c 



( 24 ) 
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one for the turbulence dissipation rate. This purely 
empirical formula is given as: 

p|£ +p y.Ve=V*(-^LVe) +c, \iA*-c t p-^ (25) 

dt Pr t * 2 T k * 2 k 

In both of these equations several similarities 
appear. The left side of both equations equate to the 
particle rate of increase of the variable. The first term on 
the right side is the diffusion rate while the second and 
third terms are the production and dissipation, respectively. 

This turbulence model proves more successful in the 
outer regions of the flow field as opposed to the inner region 
where convective transport is relatively unimportant. For 
this region wall functions have been developed for the area in 
the range of 30 < y + < 200 where y + is defined as 

i 

+ |/P„) 2 (26) 

v „ 

with T w being the shear stress, p w the density and V w 
representing the kinematic viscosity all evaluated at the 
wall . Employing these wall functions in the near wall region 
eliminates the need for a much finer grid in this area. 

B. THE GENERAL DIFFERENTIAL EQUATION 

Patankar [Ref 34] and Patankar and Spalding [Ref 35] 
recognized the relationship of the terms in the model 
transport equations and advocated a procedure of numerical 
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analysis of fluid flow which involves a general differential 
equation. Their equation is 

4 - (p<t>) +v-(pu<t>) =v-(rv<t>) +5 (27) 

dfc 

where § represents the dependent variable and T and S refer to 
the diffusion coefficient and the source term, respectively. 
These quantities are distinctive for the (J) being analyzed. 

The realization that the applicable transport equations of 
momentum, the energy, turbulence kinetic energy and 
dissipation and related coefficients can be thought of as 
particular cases of the general <)) equation is an important 
time-saving step. Now for the solution of a fluid flow and 
heat transfer problem, only the numerical solution of equation 
27 need be investigated. By constructing the applicable 
computer algorithm based on a pressure correction scheme [Ref 
34] in such a manner as to repeatedly solve this same equation 
with an appropriate definition for <)) and the pertinent values 
of r and S, considerable computational time and memory savings 
are realized. 

C. BOUNDARY CONDITIONS 

1. Boundary Conditions at a Wall 

No-slip boundary conditions are used for the 
velocities at a solid wall, where they are all set equal to 
zero. Wall temperatures are either explicitly set or the 
temperature gradient at them is specified. When the flux at 
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the wall is set equal to zero the wall assumes the properties 
of an adiabatic wall. 

2. Boundary Conditions at Inflow 

The conditions at the inflow include all the velocity, 
pressure, temperature and turbulence intensity. For this 
investigation these values were all specified at the inlet 
condition . 

3 . Boundary Conditions at Outflow 

The outflow of this problem required only the 
specification of the pressure to atmospheric conditions. All 
other conditions at the outflow were extrapolated by the 
analysis . 
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III. THE PHOENICS COMPUTATIONAL PROGRAM 



A. HOW THE PHOENICS COMPUTATIONAL PROGRAM WORKS 

Phoenics is a computer code which simulates fluid-flow, 
heat-transfer , chemical-reaction and related phenomena. It is 
capable of solving three-dimensional, unsteady and turbulent 
flow problems. Relying on the validity of the assumptions 
used to generate the input, it has the ability to simulate the 
effects of interaction between such phenomena and the 
surrounding environment. Phoenics was used throughout this 
investigation. For each dependent variable (|) , Phoenics 
generates as many difference equations as cells within the 
domain. Thus, the total number of equations Phoenics must 
solve is x times y times the number of dependent variables 
present in the domain. Not only are these equations numerous, 
but they are also very strongly coupled. In addition they are 
also characteristically non-linear due to the functional 
dependence of the coefficients and sources upon the auxiliary 
and dependent variables. To solve these non-linear equations, 
Phoenics employs an iterative predict-and-correct technique 
which reduces the imbalance between the left and right sides 
of each equation to a negligible level. 
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B . GRID GENERATION 



Phoenics can implement and use Cartesian, cylindrical 
polar or curvilinear or body-fitted coordinates. For two- 
dimensional Cartesian grids, specification of the length and 
height of the area of analysis is required. The area between 
the x and y boundaries can then be filled with as many grid 
lines as the solution requires. Spacing of these lines is 
either uniform, exponential (with continually increasing or 
decreasing spacing) , or in a "power" arrangement, which both 
increases and then decreases the spacing of the grid lines. 
These last two options allow the generation of a grid 
permitting clustering of the grid points closer together near 
the boundaries. When different regions are designated, the 
spacing can vary from region to region as necessary. In 
addition, by stipulating various regions in the x or y 
direction, irregular shapes can be analyzed by placement of 
zero-porosity obstructions which eliminate the zones from the 
flow field. 
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IV. BACKGROUND ON THE TURBULATOR DUCT PROBLEM 



A. EXPERIENCE ON CASTING PROBLEMS IN TURBINE BLADES 

As shown in Chapter I, internally routed cooling air flow 
proves quite effective in permitting higher internal 
operating temperatures for gas turbine engines particularly if 
turbulators are used in the coolant passages. The size, shape 
and orientation of the duct and the geometry of the 
turbulators within these ducts is crucial to the convective 
removal of heat from the blades. As many investigations 
showed, altering the aspect ratio of the duct or the 
turbulators, the spacing of the turbulators, their 
orientation, the turbulator pitch-to-height ratio and their 
shape, greatly impacts their effectiveness in promoting heat 
transfer . 

With the continual emphasis on reduction of weight for 
today's aircraft and their components, engine dimensions, both 
externally and internally, are also decreasing. For this 
reason the casting process by which turbine blades are 
manufactured is placing restrictions on the size and shape of 
the turbulators which can be employed. During the turbine 
blade casting process one of the problems encountered involves 
the inability of the molten materials to totally fill the 
grooves forming the turbulator. 
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A second problem is the tolerances provided by the 
manufacturing process. Current turbine sizes require 
turbulators measuring approximately 0.010" in height. The 
lower limit for associated dimensional tolerances is 0.002" - 
0.003". Thus a turbulator designed for a height of 0.010" can 
have a variance of up to 30%. This difference can 
significantly effect the heat transfer coefficient for an 
engine where a 0.010" turbulator presents a high blockage 
ratio . 

B. AN EXPERIMENTAL INVESTIGATION INTO A TURBULATED DUCT 
1. Background and Duct Construction 

Taslim and Spring [Ref 7] performed a series of 
experiments to determine the variation in heat transfer caused 
by the geometrical changes characteristic of small cast air- 
foils. Their method involved using a 46 inch long channel 
constructed of plexiglas on three sides, with the fourth side 
made of a polyurethane slab incorporating a series of heaters 
and liquid crystal sheets. The width of the heated wall was 
3", which fixed that dimension for the duct. The length of 
the side walls were changed to alter the aspect ratio during 
the experiment . Experiments were run with aspect ratios for 
the duct alternating between 0.50 and 0.55. 

The first 14 inches of the duct had no turbulators and 
simulated the unturbulated area at the root of the airfoil, 
while the final 32 inches contained a staggered array of 
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turbulators on the heated face and the opposing plexiglass 
wall. During the experiment, turbulator height , e, center-to- 
center pitch, S, and turbulator width, w, were varied. 
Blockage ratios, e/D h where D h is the hydraulic diameter, ran 
between 0.15 and 0.285, while pitch-to-height ratios. S/e, 
varied between 4.27 and 14.2. 

2 . Conduct of the Experiment 

After calibrating the liquid crystal sheets in a water 
bath so a given temperature produced a known color, and with 
the desired turbulator and duct configuration established, air 
flow was introduced into the duct with a set mass flow rate. 
This combination of duct geometry and mass flow rate 
established the Reynolds number for the test run. 

Selecting an area approximated by x/D h = 16, where x 
is the distance from the test section entrance and D h is the 
hydraulic diameter of the unturbulated section, the heat flux 
was varied by increasing or decreasing the power on the wall 
heaters until a reference color was obtained on the downstream 
side of a turbulator near the selected area. A picture of 
this area was taken and the power altered until the same 
reference color was achieved in the next adjacent area. This 
process was repeated until the entire area between the pair of 
turbulators was covered. 

The final experimental step required digitizing each 
picture and measuring the area covered by the reference color. 
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This allowed production of an area-weighted average heat 
transfer coefficient. 

3. Results 

The first result desired from the experiment was to 
determine the impact of alleviating the problems caused by the 
aforementioned non-fill difficulty by reducing the aspect 
ratio of the turbulator, AR t . Taslim and Spring [Ref 7] found 
a reduction of AR t from two to one reduced the Nusselt number 
on the turbulated wall by 33% for high blockage ratio 
turbulators (e/D h =0.22) with a spacing of S/e of 5 . As 
expected, a reattachment point formed directly downstream of 
the turbulators. This was the point of maximum heat transfer 
coefficient, which then decreased monotonically in the 
downstream direction until it approached the next turbulator, 
where it started to increase again. The conclusion from this 
analysis was that an optimum spacing between turbulators 
existed which allowed sufficient length for the flow to 
reattach, and then a minimum of spacing thereafter to the next 
turbulator. The conclusion was that, for a square turbulator 
( AR. = 1) with a relatively high blockage ratio, the optimum 
S/e spacing was approximately 8.5. 

After this determination, the effects of rounding the 
turbulators at the tops and decreasing the sharpness of the 
fillets were investigated. Testing was performed on four sets 
of turbulators. Two sets were manufactured, both with an AR t 
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of 1.0, with one set having sharp fillets and sharp corners 
and the other having rounded fillets and rounded corners. Two 
similar sets were constructed with an AR C of 2.0, with the 
sharp or rounded fillets and corners. As expected the 
rounding decreased the Nusselt numbers on the turbulated walls 
for both AR t ' s , but the AR t = 2.0 configuration produced a 17% 
reduction over the sharp cornered turbulators with the same 
AR t , while the difference for the turbulators with an AR t of 
1.0 was only 5%. Concurrently incorporated was a reduction of 
8.6% of surface area for the turbulator with an AR t of 2 and 
a 7.15% reduction for the other set. The contrast was 
attributed to the difference in spacing between the sidewalls 
of the adjacent turbulators. Although the S/e was held 
constant, the lower AR C turbulators were wider and hence less 
distance existed between sidewalls. Thus for higher aspect 
ratio turbulators, the flow had sufficient length to reattach 
after the turbulator, while for the lower aspect ratio, 
reattachment did not occur. In the latter case (tip-to-tip 
turbulator flow) there was less dependence on turbulator 
geometry . 

From these results an optimum turbulator design could 
be attempted. For manufacturability, the corners would be 
rounded and the base would be wider. For maximum heat 
transfer the tip would be kept narrow. The resulting rounded 
trapezoidal shape is termed a "jersey-barrier". At S/e of 
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five the jersey-barrier shaped turbulators exhibited the same 
Nusselt number as high AR t turbulators. 

Another area investigated was the impact of 
manufacturing tolerances on the heat transfer characteristics 
of turbulated ducts. Using a jersey-barrier shaped turbulator 
with a nominal height of 0.010" and a tolerance of ± 0.003" on 
S/e of 10, a ±40% deviation in the turbulated wall Nusselt 
number was experienced when dimensions were allowed to 
approach their minimum and maximum tolerance values. Reducing 
the S/e ratio to 5, caused a variation from the nominal case 
of +22% and -13% from the maximum to minimum profiles. 

The final area examined was the interaction between 
the blockage ratio, the pitch-to-height ratio and the 
turbulator aspect ratio. It was found that the decrease in 
Nusselt number associated with a reduction in e/D h could be 
completely offset by increasing the AR C and S/e. In addition, 
with a constant S/e close to the optimum level and a moderate 
blockage ratio, decreasing the AR C from 1.0 to 0.625 resulted 
in only a 7% drop in the Nusselt number. This was 
characteristic for geometries with the S/e closer to the 
optimum value (8.5) and a lower value for the blockage ratio. 
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V. 



PHOENICS SIMULATION OF THE TURBULATOR DUCT EXPERIMENT 



The Phoenics program was used to compute the results of 
the previous turbulated duct experiment. The analysis was 
divided into two basic steps. The first was to generate a 
rectangular straight duct with one wall heated. By 
controlling the Reynolds number, through varying inflow 
velocity, the Dittus-Boelter correlation of NU S = . 0243Re°‘ 8 Pr 0 ‘ 4 
was verified [Ref 31]. The simulation was two-dimensional. 

The second phase was the generation of a turbulated duct 
to reproduce experimentally obtained results through numerical 
analysis. Again a rectangular two-dimensional straight duct 
was used with a staggered arrangement of turbulators on the 
top and bottom wall. The bottom wall was heated. 

A. DUCT AND FLOWFIELD SETUP 
1. Grid Generation 

a. Unturbulated Square Duct 

The first duct to be modeled was 30.0m long. The 
height of the duct was 1.0m and it contained no turbulators. 
The bottom wall was heated. 

In the x direction the entire duct was considered 
one region. It was divided into 100 grid points which were 
stretched in the x direction using a factor of 1.10. This 
configuration clustered the grids near the inlet. By doing 
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so, the grid cells were square at the entrance to the region, 
which insured that the flow achieved its fully developed 
profile prior to entering the zone where it was analyzed. 

In the y direction there was only one region. 
There were 50 grid points along the y axis and they were 
clustered along the top and bottom walls. 

b. Turbulated Rectangular Duct 

In the second phase two ducts were generated. The 
first duct was 1.17m long with the initial 0.356m 
unturbulated . The turbulators had an aspect ratio, AR t , of 
1.0, with a height, e, of 0.0127m and a turbulator center-to- 
center pitch, S, of 0.0635m. All the turbulators were square, 
with sharp corners and fillets. The first turbulator was on 
the top wall and all were positioned 90° to the flow. The 
second duct was 0.654m long, with the same unturbulated region 
at the inlet. The duct aspect ratio remained the same. Ten 
turbulators of the same AR t , S, and e were again positioned 
90° to the flow, with the first one on the top wall. The 
bottom wall of both ducts in the second phase was heated. 

The longer duct was divided into 30 regions along 
the x axis. It included 175 grid points which were stretched 
by a factor of 1.2. Region two and every fourth region there- 
after (6,10,14, . . .) through number 2 6 were the regions designa- 
ted for the turbulators on the top wall. Similarly region 
four and every fourth region following through number 28 were 
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the regions demarcating the locations of the bottom 
turbulators . All of the regions designating the turbulators 
were 0.0127m long and contained seven grid points. The odd 
numbered regions through number 29 specified the spaces 
between the upper and lower turbulators . Each was composed of 
nine grid points and was 0.0191m in length. The final region 
was 0.3 69m long and modeled outflow from the duct. No 
turbulators were located here and it contained 101 grid 
points . 

In the vertical direction three regions were 
specified. The first region was from the bottom wall up 
0.0127m, which corresponded to the height of the turbulators. 
It contained 15 grid points. The second region extended from 
the top of the first up 0.0165m. Twenty-five grid points were 
allocated in this region. The final region in the y direction 
simulated the top wall turbulators. Like region one it was 
0.0127m tall and held 15 grid points. Figure 1 shows a 
schematic of the longer duct with an enlarged segment which 
depicts the grid at the end of the first region and the first 
two turbulators in the turbulated region. 

This grid most closely resembled the experiment, 
but did not give adequate grid resolution. An additional duct 
was generated, which had a finer mesh size. This duct was a 
shortened version of the fully turbulated one. As will be 
explained, investigation revealed that the shorter length did 
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not significantly alter the flowfield prediction and did allow 
a greater concentration of grid points throughout the duct. 
This final duct also included an unturbulated inlet region, 
whose dimensions and grid points were the same as the previous 
duct. The shortened turbulated region measured 0.298m long 
and consisted of ten turbulators arranged in a staggered array 
again with the first one on the top wall. Each turbulator has 
fifteen grid points from front to back and there are 20 grid 
points in the space between each turbulator. This duct was 
investigated using a uniformly spaced grid and stretching 
factors of 1.05, 1.10 and 1.20 all of which clustered the 

grids near the front and back of the turbulators. There were 
325 grid points in this region and a total of 500 in the 
entire horizontal or streamwise direction. 

The vertical direction was constructed the same 
way as the previous duct with three regions of fifteen, 
twenty-five and fifteen grids each. Grid spacing here also 
employed both stretched and uniformly spaced points. Figure 
2 shows the same area for this duct as in Figure 1 . Figure 2 
uses the 1.2 level of grid spacing on the increased level of 
grid points. 

2 . Boundary Conditions 

Analyses were run at a number of different Reynolds 
numbers and turbulence levels. The Reynolds number is defined 
as 
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Figure 1 Schematic of Long Duct With Enlargement of the 
First Two Turbulators Placed in a Uniformly Spaced Grid 
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( 28 ) 



where D h is the hydraulic diameter. For a rectangular duct, 
D h is defined as 4A/P where A is the area of the duct and P is 
the perimeter. For a two dimensional analysis a length scale 
is used. For the purposes of this investigation the D h used 
was the value for the actual D h of the duct as configured in 
Taslim's and Spring's experiment. Using the values for 
density of 1.161 kg/m 3 and for viscosity of 1.589 X 10' 5 kg/m- 



s, Table 1 gives the inlet velocities in the horizontal 
direction associated with the Reynolds numbers which were 
investigated . 



Reynolds number of 10000. The turbulent kinetic energy is 
calculated as follows 



where Tu is the turbulence intensity and u is the inlet 
velocity, the values for the appropriate turbulence levels to 
keep k approximately equal at the inlet for all Reynolds 
numbers is also given in Table 1. 



A value of 2 



% 



inlet turbulence was specified for a 



k= ( Atu) x [ u] 2 



( 29 ) 
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REYNOLDS NUMBERS, VELOCITIES AND TURBULENCE LEVELS 



INVESTIGATED 



Reynolds 


Velocity 


Turbulence 


Number 


(m/s) 


(%) 


10000 


2.531 


2.000 


15000 


3.796 


0.889 


20000 


5 . 062 


0 . 500 



Table 1 

The inlet flow temperature was held at 300° K for all 
runs. The bottom wall was heated with a constant temperature 
of 400° K, and the top wall was specified as adiabatic. At 
the outlet the external static pressure was specified. 

The turbulators were obstructions in the flow and so 
they were specified as non-conducting, adiabatic obstructions. 
In addition they were given the characteristics of possessing 
friction and having zero-porosity. No gravity forces, heat 
sources or internal plates were simulated in this analysis. 
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VI . COMPUTATIONAL RESULTS 



A. UNTURBULATED PLANAR DUCT 

As an initial test case an unturbulated square duct was 
evaluated by the Phoenics computational program. For fully 
developed laminar flow, a parabolic velocity profile develops 
with a maximum velocity, u roax =2u mean . For turbulent flow the 
velocity profile is much flatter in the center due to the 
mixing motion and momentum exchange and u raax =l . 2u raean . 

The values of The Reynolds numbers for this investigation 
were 10000 and 20000. The geometry of the chosen duct 
provided a 30:1 length to height ratio. This allowed a 
sufficient length for the flow to develop fully prior to being 
investigated. The zone chosen for the analysis of the flow 
was located between 20 and 25 diameters downstream of the 
inlet . 

Figure 3 is a profile of the flow located 25 diameters 
from the inlet. As can be seen it presents the classic shape 
associated with turbulent fluid flow in a duct. The maximum 
value of the velocity of the flow occurs at mid-duct and is a 
value of 0.1181 m/s. Using the Phoenics results file the mean 
average for velocity at this point in the flow was 0.998 m/s. 
The ratio is approximately 1.18, which is indicative of 
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Figure 3 - Flow Profile 25 Diameters Downstream in a Square 
Unturbulated Duct 



turbulent flow. A plot of the skin friction coefficient, c £ , 
where 







( 30 ) 



T w is the shear stress along the wall and S is the area, on 
the bottom wall (Figure 4) shows the increase associated with 
the viscous effects as the flow progresses along the duct. 
The value for this shear stress along the top and bottom wall 
for the area of fully developed flow is 8.175 X 10‘ 5 kPa . 
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B. TURBULATED RECTANGULAR DUCT 
1. Full Turbulated Duct 

A solution for the turbulated full duct was generated 
first. As stated previously this analysis was run with a 
relatively coarse grid in order to simulate Taslim's [Ref 7] 
experiment. This run included a Reynolds number of 10000 and 
converged at approximately 1900 iterations. 

As shown in Figure 5, at both the first and second 
turbulators the flow generated a stagnation point on the 
upstream face of the turbulator and recirculation occurred 
ahead of the turbulator. Figure 6 is a plot of the velocity 
vectors for the area of the duct adjacent to the second 
turbulator. The flow did not recirculate upstream of 
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Figure 5 - x-Direction Velocity Vectors at Turbulators 1 and 
2 
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Figure 6 - Velocity Vectors Between Turbulators 2 and 4 

Using a Constant Grid Spacing 
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turbulators three and four, because of the influence of the 
upstream turbulators two and three respectively. This flow 
pattern became repeatable from this location onward until the 
end of the turbulated section. 

The position of the reattachment point downstream of 
the turbulator was dependent on the location within the duct. 
As the flow progressed down the duct, the reattachment point 
moved closer to the downstream face of the turbulator 
eventually reaching a fixed distance downstream of the 
turbulator. The area between turbulators eight and ten was 
the first to approximate this "periodic" flowfield. Figure 7 
shows the vectors as they appeared for the flow in this area 
and Figure 8 plots the x-direction velocity components for the 
first grid point away from the lower surface. 

Another area of interest was the behavior of the 
freestream kinetic energy (k) of the turbulence. The mid-duct 
value for turbulent kinetic energy was plotted for the entire 
length. Figure 9 shows an initial free stream decay of 
turbulence in the unturbulated section, after which the 
turbulence increased dramatically at the first turbulator and 
experienced additional increases at each subsequent one. 
Initially each gain was significant, but eventually they 
became smaller in magnitude. After the eighth turbulator a 
near constant value was approximated for the free-stream 
turbulence. This asymptotic value was two orders of magnitude 
over the minimum values experienced in the unturbulated 
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Figure 7 - Velocity Vectors Between Turbulators 8 and 10 

Using a Constant Grid Spacing 
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Figure 8 - x-Direction Velocity Component Between 

Turbulators 8 and 10 Using a Constant Grid Spacing 



section. After the last turbulator, the flow entered the 
unturbulated outlet portion of the duct and turbulence once 
again decreased due to viscous dissipation. 




Figure 9 - Turbulent Kinetic Energy for a Full Turbulatec 
Duct 
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Using the above data, the tenth turbulator was 
employed as the cutoff point for the full investigation. 

2 . Abbreviated Turbulated Duct 

a. Uniformly Spaced Grid Simulation 

Using a uniformly spaced grid and a duct shortened 
to ten turbulators, a simulation was run at a Reynolds number 
of 10000 with a turbulence intensity of 2.0%. The development 
of the flow was consistent with the full length duct. Both 
the first turbulator (on the top wall) and the second one (on 
the bottom wall) developed regions of recirculation prior to 
them. Figure 10 is a plot of the velocity vectors for the 
area prior to turbulator two and it shows the flow pattern 
with a distinct stagnation point. 

The flow for the area downstream of turbulator 
two, however was markedly different from the structure found 
using the full length duct. Here the entire area between the 
second and fourth turbulators exhibited recirculation of the 
flow. No reattachment point along the bottom wall existed. 
Figure 11 is a plot of the vectors in this area and shows the 
total reverse flow. A plot of the horizontal velocity at the 
first grid line off the bottom wall for this region shows the 
negative value of this component in Figure 12. 

After the fourth turbulator, the flow along the 
bottom wall begins to exhibit the characteristics of the pre- 
vious duct. The flow over the turbulator turns down to the 
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Figure 10 - x-Direction Velocity Vectors Prior 

Turbulators 1 and 2 Using a Constant Grid Spacing 



to 
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Figure 11 - Velocity Vectors Between Turbulators 2 and 4 
Using a Constant Grid Spacing 
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reattachment point establishing a recirculation zone 
downstream of the turbulator. The flow then turns toward the 
next downstream turbulator. On the upstream face of this 
turbulator no recirculation occurs and the flow proceeds over 
the top with no separation on the upper surface. 
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Figure 12 - Horizontal Velocity Component Between 

Turbulators 2 and 4 Using a Constant Grid Spacing 



Using data generated in the Phoenics results 
files, the area between the eighth and tenth turbulators was 
investigated further. Figures 13 and 14 show the horizontal 
and vertical velocity components respectively of the flow in 
this area along the first grid line off the bottom wall. The 
horizontal velocity shows a minimum value of -0.51 m/s in the 
recirculation zone with the reattachment point located 0.014m 
downstream of the turbulator. The flow accelerates for 0.027m 
where it reaches a maximum of 1.9 m/s. After the reattachment 



54 




point all values for the flow are positive, verifying the lack 
of a recirculation region on the face of the tenth turbulator. 

Figure 14 illustrates the vertical velocity 
component's maximum values are located in those areas immedi- 
ately after and immediately prior to the turbulator faces. 




Fignre 13 - Horizontal Velocity Component Between 

Turbulators 8 and 10 Using a Constant Grid Spacing 



Figure 15 shows a plot of the calculated 'Stanton number versus 
location in the duct. The first peak is located adjacent to 
the downstream face of turbulator eight and coincides with the 
maximum vertical component of flow and a near zero value of 
horizontal velocity. Due to the location of this peak 
adjacent to the wall of the turbulator, predictions of global 
heat transfer were not consistent with those found in the 
experiment. This will be examined closer in the next section 
using grid stretching to cluster more grids near the wall. 
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8 and 10 Using a Constant Grid Spacing 




Figure 15 - Stanton Number Between Turbulators 8 and 10 
Using a Constant Grid Spacing 

The next three peak Stanton numbers occur between 
0.010m and 0.011m from the downstream face of the turbulator. 
When compared to the velocity component plots, two distinct 
peaks are found to straddle the reattachment point with one 
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rather elongated peak occurring over the reattachment zone. 
The final peak occurs near the upstream face of the tenth 
turbulator. This maximum is located at the position of the 
maximum vertical component of the flow and where the 
horizontal component is nearing zero. 

The effect of the fluid flow and heat transfer 
within the duct is shown in Figure 16 which .is a plot of the 
wall temperature between the two turbulators . The maximum 
temperature develops near the downstream face of the eighth 
turbulator. The temperature first decreases and then in- 
creases to a relative maximum prior to the reattachment point. 




Figure 16 - Temperature Between Turbulators 8 and 10 Using 
a Constant Grid Spacing 



At the reattachment point a relative minimum temperature re- 
sults from the predominantly downward flow. Another relative 
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maximum occurs on the downstream side of the reattachment 
point followed by the minimum temperature region. This region 
corresponds to the area where the horizontal velocity 
component is, in general, one order of magnitude greater than 
the vertical component . A relative temperature maximum arises 
near the upstream face of the tenth turbulator, corresponding 
to the relative Stanton number maximum in the same location. 

The development of the kinetic energy of 
turbulence is plotted in Figure 17. An initial decrease in 
the unturbulated area of the duct approaches a minimum of 4.3 
X 10' : '. As the flow enters the turbulated region the effects 
of each turbulator is observed. The characteristic increases 
experienced by KE at each turbulator initially are on the 
order of 3.5 to 4.0. Eventually they decline to approximately 
0.1. The total increase of KE approximates two orders of 
magnitude . 

b. Simulation Using Grid Stretching 

Using the same duct as above but with a 1.2 factor 
to stretch the grid in all regions in both the x and y 
direction, another simulation was run. Again the flow was run 
at a Reynolds number of 10000 with an inlet turbulence 
intensity of 2.0%. The initial development of the flow 
remained the same. The characteristic recirculation regions 
formed prior to both turbulators one and two. The area 
between turbulators two and four again contained totally 
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Figure 17 - Turbulent Kinetic Energy for the Shortened Duct 

separated flow with no reattachment point. Subsequent flow 
developed in a pattern similar to the flow associated with a 
uniformly spaced grid. A reattachment point developed after 
the turbulator with an area of recirculation upstream of it 
between the turbulator and the reattachment point. No 
recirculation on the upstream face occurred for any flow after 
the second turbulator. 

Again using the Phoenics results file the area 
between turbulators eight and ten can be examined closely. 
Figures 18 and 19 show the horizontal and vertical components 
of the flow velocity respectively in the analyzed regions. 
Figure 18 indicates this grid generated a larger recirculation 
zone than the uniformly spaced grid. The stretched grid 
developed a recirculation zone of 0.017m versus the uniform 
grid's 0.014m. The stretched grid also generated smaller 
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Figure 18 - Horizontal Velocity Component Between 

Turbulators 8 and 10 Using a 1.20 Stretching Factor 




Figure 19 - Vertical Velocity Component Between Turbulators 
8 and 10 Using a 1.20 Stretching Factor 

maximum and minimum values for the vertical component of the 
velocity although the locations were relatively similar. 
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The plot of the Stanton number between the two 
turbulators (Figure 20) showed considerable differences when 
coupled with the one for the uniform grid (Figure 15) . 




Figure 20 - Stanton Number Between Turbulators 8 and 10 
Using a 1.20 Stretching Factor 

The stretched grid solution located a relative 
maximum adjacent to the downstream face where the vertical 
component of the velocity was a maximum, and the horizontal 
component was near zero. It also located two relative 
maximums on either side of the reattachment point as 
determined from Figure 16. The next point identified is the 
maximum for the Stanton number in the area between these two 
turbulators. It is positioned just downstream from the 
minimum vertical velocity component, approximately 0.021m from 
the downstream face of turbulator 8. Its effect is 
discernable from Figure 21, the plot of the temperature 
distribution between the two turbulators. 
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Figure 21 - Temperature Between Turbulators 8 and 9 Using a 
1.20 Stretching Factor 

The maximum temperature is located precisely where 
the maximum heat transfer occurs, as indicated by the Stanton 
number. The final relative . peak for the Stanton number is 
located proximate to the upstream face of the tenth turbu- 
lator. As in the previous solution it is collocated with the 
maximum value for the flow's vertical component of velocity 
and where the horizontal component is approaching zero. 

Since the maximum Stanton number and temperature 
occurs near the position on the grid where the clustering of 
the grid points for the upper turbulator is greatest, it ap- 
pears that this maximum is grid related, rather than a valid 
depiction of the actual situation. Another variation of this 
configuration was run using a grid stretching factor of 1.10. 
As seen in Figure 22, the peak appearing in Figure 20 is not 
present for this new configuration. The new peak is located 
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Figure 22 - Stanton Number Between Turbulators 8 and 10 

Using a 1.10 Stretching Factor 

on the upstream side of the attachment point and is half the 
magnitude of that experienced for a stretching of 1.20. 

A final attempt to examine this grid generated 
anomaly was made using a grid stretching factor of 1.05. 
While this tended to cluster the grid closer to the downstream 
face of the turbulator, which appeared to rectify the spike 
from the uniformly spaced grid, it did not create an artifi- 
cial spike in the area adjacent to or in the reattachment 
zone. Figure 23 is the Stanton number plot for the flow be- 
tween turbulators eight and ten for this grid configuration. 
The most obvious difference between this and the previous 
cases is the reduction in the magnitudes of the peaks. All 
the peaks in this plot which are located at the same approxi- 
mate location as the peaks associated with the previous grids 
are lower by an order of magnitude. The peak adjacent to the 
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downstream face of the turbulator is also approximately an 
order of magnitude smaller than the peak associated with the 
uniformly spaced grid. This analysis confirms that the 
solution for this area is highly grid dependent. 
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Figure 23 - Stanton Number Between Turbulators 8 and 10 

Using a 1.05 Stretching Factor 



The plot of the temperature for the area between 
turbulators eight and ten using a 1.05 stretching factor 
(Figure 24) shows that the maximum temperature is associated 
with the relative maximum Stanton number adjacent to the 
downstream face of the turbulator. The other three relative 
maximums for temperature are also associated with relative 
maximums of the Stanton number. Contrary to all previous 
investigations a relative minimum of the temperature occurs at 
a relative maximum Stanton number. 

Defining shear stress according to equation 13 and 
plotting it for this portion of the duct results in Figure 25. 
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Figure 24 - Temperature Between Turbulators 8 and 10 Using 
a 1.05 Stretching Factor 



Here the plot shows that the value of the shear stress at the 
point of the relative maximum Stanton number is zero. 





Figure 25 - Shear Stress Between Turbulators 8 and 10 Using 
a 1.05 Stretching Factor 



This duct's turbulent kinetic energy plot (Figure 
26) is very similar to that for the previous duct. Original 
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values are the same with the same initial decrease 
experienced. Although they both reached the same level of 
turbulence, the stretched grid appears to take longer to 
develop this value. The final increases, however, are still 
more than an order of magnitude smaller than the initial ones 
and can be considered to be an asymptotic condition. 




Figure 26 - Turbulent Kinetic Energy for the Abbreviated 
Duct Using a 1.05 Stretching Factor 

Another investigation which was done was the 
determination of the effects of varying the turbulence inten- 
sity at the inlet. Figures 27 and 28 show plots of the KE at 
a mid-duct location for a Reynolds number of 20000 with the 
initial turbulence intensity being 2.0% and 0.5% respectively. 
As expected the starting and development areas of the duct for 
both conditions are different. However once the flow enters 
into the turbulated area, the flow patterns and the levels of 
turbulence start to converge. From this it can be concluded 
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that the level of turbulent kinetic energy at an asymptotic 
location in a duct is unaffected by the level of turbulence 
intensity far upstream. 





Figure 28 - Turbulent Kinetic Energy for a Reynolds Number 
of 20000 and a Turbulence Intensity of 0.5% 
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VI. COMPARISON OF RESULTS WITH EXPERIMENT 



A. UNTURBULATED SQUARE DUCT 

The results of the unturbulated square duct can be 
compared with the results of the Dittus-Boelter equation which 
states 

Nu=0 . 0243Pr°- 4 Re°- 8 ( 31 ) 

where Nu is the Nusselt number. This equation is valid when 
the Prandtl number is between 0.7 and 160.0 and the Reynolds 
number is greater than 10 4 . From Phoenics, the steady state 
Stanton number was calculated for Reynolds numbers of 10000 
and 20000. This is a measure of the heat transfer occurring 
at the bottom wall of the duct. This number was then 
converted to the Nusselt number. Figure 29 is a plot of the 
Dittus-Boelter equation (dashed line) and the experimental 
data between Reynolds numbers of 10000 and 20000. 

The disagreement between the Phoenics generated values for 
Nusselt number and the Dittus-Boelter correlation is 
considerable. The slope is opposite between the two and the 
magnitude is one full order higher for the Phoenics generated 
data. Operating at the lower end of the effective region for 
the Dittus-Boelter correlation and grid resolution problems 
are possible causes of this discrepancy. 
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Duct 

B. TURBULATED RECTANGULAR DUCT 
1. Constant Grid Spacing 

The plot of the Stanton number for the uniformly 
spaced grid at a Reynolds number of 10000 (Figure 15) is 
dominated by the peak near the downstream wall. This is 
specific to this test case and does not appear on the plots 
for Reynolds numbers of 15000 or 20000. This adversely 
affects the entire slope of the numerically generated Nusselt 
number solutions . 

Since each grid is a constant width, the average 
Stanton number over the entire zone can be resolved by deter- 
mining the value for the Stanton number at each grid point and 
averaging it. Plotted on log-log paper with the results of 
Taslim and Spring [Ref 7], the results of the uniformly spaced 
grid show a negative slope, while that of the experimentally 
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generated data is positive. In addition, the experimental 
data corresponds closely to the slope of the Dittus-Boelter 
correlation line. A plot of these two lines is shown in 
Figure 30 with the experimental data being depicted with a 
dashed line. 




Figure 30 - Nusselt Number vs Reynolds Number for a 

Uniformly Spaced Grid 

2. Grid Spacing Using Stretching 

Three different factors of stretching were used 
throughout this analysis: 1.20, 1.10 and 1.05. Since the 

grids are variably spaced, the average Stanton number requires 
each individual Stanton number be multiplied by the width of 
its grid. All of these values must then be added together and 
then divided by the total width to determine the average. 
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The plots of the Stanton numbers of the 1.20 and 1.10 
stretched grids are dominated by exceeding high artificial 
peaks (Figures 20 and 22) . When converted to Nusselt numbers, 
averaged with the other values and plotted against Reynolds 
numbers, these artificially high values induce negative slopes 
on the resulting curves. This runs contrary to not only the 
results of Taslim, but also the Dittus -Boelter correlation and 
the Colburn correlation which states Nu=f /2Pr 1/3 Re where f/2 is 
a dimensionless shear stress term [Ref 36] . All of these 
clearly indicate an increase of Reynolds number should cause 
an increase of Nusselt number. A plot of the resultant graph 
for the stretching of the grid by factors of 1.10 and 1.20 are 
shown in Figures 31 and 32 respectively. 




Figure 31 - Nusselt Number vs Reynolds Number Using a Grid 
Stretching Factor of 1.10 
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REYNOLDS NUMBER 



Figure 32 - Nusselt Number vs Reynolds Number Using a Grid 
Stretching Factor of 1.20 



The stretching factor of 1.05 produced no such 
anomalies. The values for Stanton number were more 
consistent, with no one peak exceeding the others by an order 
of magnitude or more. The plot of this resultant curve is not 
similar to the slope of that generated by Taslim and Spring. 
As seen in Figure 33, with its minimal slope, it more closely 
approximates a constant Stanton number instead of an 
increasing one with increasing Reynolds number. 

The calculation of the Stanton number by the Phoenics 
algorithm is a blind operation which does not reveal the 
equations or values used. To determine if the variance in the 
Phoenics generated Nusselt numbers are the result of a grid 
resolution problem, the Nusselt number was hand calculated 
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Figure 33 - Nusselt Number vs Reynolds Number Using a Grid 
Stretching Value of 1.05 



from the values of component variables output by Phoenics and 
compared with the values computed in the code. 

First the heat flux, q 0 '', is defined as 



Qq = (*>P 



C»C p k 2 



dT 

dy 



(32) 



where K is the thermal conductivity of air at 300°K, is a 
constant in the k-8 turbulence model equal to 0.09, C p is the 
specific heat of air at constant pressure and a t is the 
turbulent Prandtl or Schmidt number [Ref 37] . 

Once the value of the heat flux is determined, a local 
heat-transfer coefficient, h x , is specified as 
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( 33 ) 




where t 0 is the specified temperature at the bottom wall. The 
t„ temperature is located at mid-duct. Finally the local 
Nusselt number can be generated from 

h v x , _ . . 

Nu=—Z— ( 34 ) 

K 

Here the x value is taken as the width of the cell for which 
the temperature, kinetic energy and dissipation of turbulence 
are being considered . [Ref 38] 

A plot of the calculated Nusselt number versus the 
distance from the inlet is shown in Figure 34. The dashed 
line plot is the Phoenics generated values. A very striking 
resemblance exists between the two graphs. All but one of the 
maximums from the Phoenics results are identifiable for the 
calculated values . The one value which does not appear on the 
chart of the calculated values is the spike at the reattach- 
ment point. The indication is this spike is grid generated. 

As a final comparison, the values for the average 
calculated Nusselt numbers for Reynolds numbers of 10000 and 
20000 were plotted versus these Reynolds numbers in Figure 35. 
Also plotted were the values obtained from Phoenics for a 1.05 
stretching factor solution (solid line) and the experimental 
values (dashed and dotted line) . The slope for the calculated 
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DISTANCE FROM INLET (m) 



Figure 34 - Calculated and Phoenics Generated Values for 

Nusselt Number Between Turbulators 8 and 10 

values is much closer to the experiment than that produced by 
Phoenics, but the magnitude is much greater. 

Because of the close correlation between the Phoenics 
generated program and the hand calculated one, the failure to 
produce Taslim's results over a range of Reynolds numbers 
occurs due to an inability to adequately resolve the flowfield 
with this grid size. 
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Figure 35 - Nusselt Number vs Reynolds Number Using 

Calculated, Phoenics Generated or Experimental Values 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The analysis of the data presented in this investigation 
leads to a number of conclusions about airflow in a 
rectangular duct with square turbulators placed in a staggered 
array. Using a S/e ratio of 5, an e/D h of 0.235, an AR t of 1.0 
(Taslim's value of D h ) , with an e of 0.0127m, the following 
can be concluded: 



• Flow approaching the first set of turbulators will 
experience recirculation on the upstream side of only the 
first pair of turbulators. All turbulators thereafter 
will direct the flow up and over the turbulator and will 
not generate a separation point for the recirculation of 
the oncoming flow. 

• Flow over the top of the turbulators will remain attached 
and create no recirculation zones. This was true for a 
range of Reynolds numbers between 7000 and 30000 with 
turbulence intensity varying between 0.22 and 4.1. 

• No reattachment point will develop for flow over the first 

or second turbulator. The entire region between 

turbulators two and four will be recirculating flow. This 
is the only area of the duct which will experience this 
flow structure. 

• The pattern of flow between all sets of turbulators after 
the initial two on either wall will maintain the same 
characteristic structure. The flow will develop a 
recirculation region directly after the turbulator 
followed by a reattachment point. An area of acceleration 
of flow velocity toward the downstream turbulator will 
ensue with the flow traveling over the turbulator. Each 
succeeding section results in a smaller recirculation zone 
until an asymptotic dimension is reached. 
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• Maximum heat transfer occurs at the reattachment point of 

the downstream recirculation zone. In general four 
relative maximum heat transfer locations occur: on the 

downstream face of the turbulator, on either side of the 
reattachment point and on the upstream face of the 
following turbulator. 

• The kinetic energy of turbulence throughout a turbulated 
duct maintains a distinct trace. Each turbulator 
increases the level as the medium passes over it. Initial 
turbulators provide a significantly greater level of 
increase than those located further along the duct. After 
eight turbulators, an asymptotic level of kinetic energy 
is reached. 

• The asymptotic level of turbulence in a duct is 
independent of the turbulence intensity at the inlet. In 
this duct, the individual variation of kinetic energy 
caused by each turbulator increased or decreased with 
respect to a similar change in turbulence intensity. Thus 
altering the turbulence intensity did not change the 
results of the simulation. 

• The detailed results are highly grid dependent. The 
choice of uniform or clustered grids not only effects the 
results, but the power of the clustering can also impact 
them. A converged solution can exhibit grid-generated 
phenomena which have no basis in actuality and for which 
the investigator must constantly inspect. 

• This investigation was a 2-dimensional analysis of an 

experiment performed in three dimensions . The 

experimental results showed a variation in the spanwise 
direction of the duct. The effects of the walls on the 
flow caused lower values for the heat transfer coefficient 
near the sidewalls with an increase closer to the 
centerline. A 2 -dimensional analysis would not have the 
effects of the sidewalls included and would produce this 
centerline higher value. This possibly explains the 
elevated values seen in the computed results. 

B . RECOMMENDATIONS 

Phoenics is a powerful program with the ability to solve 
highly complex flow problems. Its strength lies in the ease 
with which parameter can be altered and reevaluated for a 
given problem. The present analysis was done using strictly 
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a Cartesian coordinate system in two dimensions. With 
Phoenics ' versatility an examination of alternative turbulator 
shapes and geometries is possible. Even the evaluation of the 
"jersey-barrier" shape in three dimensions is possible using 
body-fitted coordinates. 

The non-linearity of the flow and the large number of 
equations which are solved required thousands of iterations 
for a converged solution. This can take an exceptionally long 
time even when operating on a 486-processor based computer. 
A 3000 iteration solution using the configuration of this 
report took 50 hours to complete. 

Phoenics also requires the user to accept certain 
parameters and algorithms without the option to alter them. 
Although numerous areas are accessible, many of the turbulence 
parameters, wall functions and other built-in subroutines are 
not . This precludes making changes to fit the problem being 
examined . 

A full analysis of a duct with any style turbulator should 
include a three dimensional examination of the flow. The 500 
X 55 grid used here should be considered to be a minimum in 
two dimensions. The spanwise direction would add 
significantly more points in the grid and ultimately more 
equations to be solved. This will be excessively time 
consuming and may or may not provide the accuracy required. 

For these reasons the recommendation of this report is to 
use Phoenics only for the first approximations of the solution 
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to the duct problem using coarse grids. The ease with which 
programs can be entered, changed and run makes it an excellent 
candidate for this type of initial investigation. It can 
locate areas of more increased flow activity and aid in the 
delineation of zones requiring clustered gridding. Several 
iterations with increasingly finer and more clustered grids 
will assist in either supporting or refuting the initial grid 
refinements . 

Final analysis should be carried out using a program 
possessing at least as much power as Phoenics but with more 
flexibility and a faster running environment. Two such 
programs are Proteus 3-D and Overflow. 
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